2-2

CompConj[A_] := A /.Complex[0,n_]->-Complex[0,n]

Quantum Motion of a Particle in a Box

Nicholas Wheeler
Reed College Physics Department
20 August 2008

Introduction

In Part 2 we studied solutions of

-((ℏ^2)/(2m))((∂_(x,x))+(∂_(y,y)))ψ[x,y,t]=i ℏ (∂_t) ψ[x,y,t]

which vanish on the boundaries of the unit box: 0 < {x, y} < 1. We were, in particular, led to the normalized eigenfunctions

Ψ[x_,y_,m_,n_]:=2 Sin[m π x]Sin[n π y]

{{Plot3D[Ψ[x,y,1,1], {x,0,1},{y,0,1}, Ticks->False],
Plot3D[Ψ[x,y,1,2], {x,0,1},{y,0,1}, Ticks->False],Plot3D[Ψ[x,y,1,3], {x,0,1},{y,0,1}, Ticks->False]},
{Plot3D[Ψ[x,y,2,1], {x,0,1},{y,0,1}, Ticks->False],
Plot3D[Ψ[x,y,2,2], {x,0,1},{y,0,1}, Ticks->False],
Plot3D[Ψ[x,y,2,3], {x,0,1},{y,0,1}, Ticks->False]},
{Plot3D[Ψ[x,y,3,1], {x,0,1},{y,0,1}, Ticks->False],
Plot3D[Ψ[x,y,3,2], {x,0,1},{y,0,1}, Ticks->False],
Plot3D[Ψ[x,y,3,3], {x,0,1},{y,0,1}, Ticks->False]}}//MatrixForm

These we assembled into normalized functions

Ψ[x,y]=(∑_m)(∑_n)(Subscript[C, mn])Ψ[x,y,m,n]

(∑_m)(∑_n)  (Subscript[C, mn])(Subscript[(OverBar[C]), mn])=1

which we launched into motion by introduction of appropriate exponential "buzz factors"

Ψ[x,y,t]=(∑_m)(∑_n)(Subscript[C, mn])Ψ[x,y,m,n](e^(-i ω[m,n]t))

ω[m,n]=(π^2)((m^2)+(n^2))

We looked first to cases of the form

Ψ[x,y]=a[x]b[y]=(∑_m)(∑_n)(Subscript[A, m])(Subscript[B, n])Ψ[x,y,m,n]

(∑_m)  (Subscript[A, m])(Subscript[(OverBar[A]), m])=(∑_n)  (Subscript[B, n])(Subscript[(OverBar[B]), n])1

which give rise to probability densities of the factored form

P[x,y,t]=f[x,t]*g[y,t]

In such cases x and y are at all times statistically independent. We identified cases in which the expected position

{(Subscript[⟨x⟩, t]),(Subscript[⟨y⟩, t])}=(∫_0^1)(∫_0^1){x,y}P[x,y,t]dxdy

traces a Lissajous figure. And other cases in which P[x, y, t] itself has the form (at least briefly) of a localized pimple tracing a Lissagous curve.

Then, following Chen et al, we looked to cases in which  P[x, y, 0] does not factor (x and y are no longer statistically independent), but which derive interest from the fact that P[x, y, 0] assumes significant value only on a specified Lissajous curve.  It is about some properties of cases of the latter sort that I propose now to speak.

Simplest non-factorable boxed wavefunctions

The simplest way to linearly combine the moving normalized eigenfunctions

Ψ[x_,y_,m_,n_,t_]:=Ψ[x,y,m,n]*(E^(-I ω[m,n] t))
ω[m_,n_]:=(π^2)((m^2)+(n^2))

is pairwise, the general (already normalized) form being

Ψ[x_,y_,m_,n_,p_,q_,α_,β_,t_]:=Cos[α]Ψ[x,y,m,n,t]+(E^(I β))Sin[α]Ψ[x,y,p,q,t]

where α is a "mixing angle" and β is a "relative phase angle."

Ψ[x,y,m,n,p,q,α,β,t]

Ψstar[x_,y_,m_,n_,p_,q_,α_,β_,t_]:=CompConj[Ψ[x,y,m,n,p,q,α,β,t]]

Ψstar[x,y,m,n,p,q,α,β,t]

P[x_,y_,m_,n_,p_,q_,α_,β_,t_]:=Ψstar[x,y,m,n,p,q,α,β,t]*Ψ[x,y,m,n,p,q,α,β,t]

Proceeding

P[x,y,m,n,p,q,α,β,t]//Expand

ExpToTrig[4 (E^(-I ((m^2)+(n^2)) (π^2) t+I (π^2) ((p^2)+(q^2)) t-I β))+4 (E^(I ((m^2)+(n^2)) (π^2) t-I (π^2) ((p^2)+(q^2)) t+I β))]

we see that P[x, y, m, n, p, q, α, β, t] moves (oscillates) except when

(m^2)+(n^2)=(p^2)+(q^2)

The following source

Hyperlink["Integers as sums of two squares","http://www.math.mtu.edu/mathlab/COURSES/holt/dnt/represent3.html"]

indicates that

(1^2)+(7^2)=(5^2)+(5^5)=50

provides the smallest nontrivial instance of such a coincidence. But we will—in what immediately follows—have as much or more to do with "trivial" instances (of the form (a^2)+(b^2)=(b^2)+(a^5)=N) than with non-trivial ones.

Illustrative case 1:

Ψ[x,y,2,1,1,2,α,β,t]/.{α->ArcCos[(2/3)],β->ArcTan[2]}

(4/3) (E^(-5 I (π^2) t)) Sin[2 π x] Sin[π y]+(2/3)  (E^(-5 I (π^2) t+I ArcTan[2])) Sin[π x] Sin[2 π y]

(E^(I ArcTan[2]))//ComplexExpand

Ψ1= (E^(-5 I (π^2) t))( (4/3)Sin[2 π x] Sin[π y]+(2/3) (Sqrt[5])((1+2 I)/(Sqrt[5]))  Sin[π x] Sin[2 π y])

It is because (trivially)

(2^2)+(1^1)=(1^2)+(2^2)=5

that only a single time factor enters (multiplicatively) into the description of Ψ1, with the consequence that the associated probability density is time-independent:

Ψ1star= (E^(+5 I (π^2) t))( (4/3)Sin[2 π x] Sin[π y]+(2/3) (Sqrt[5])((1-2 I)/(Sqrt[5]))  Sin[π x] Sin[2 π y])

Ψ1star*Ψ1//FullSimplify

P1[x_,y_]:=(16/9) (4 ((Cos[π x])^2)+4 Cos[π x] Cos[π y]+5 ((Cos[π y])^2)) ((Sin[π x])^2) ((Sin[π y])^2)

Plot3D[P1[x,y],{x,0,1},{y,0,1}, Ticks->False]

Illustrative Case 2:

This is but a slight variant of the preceding case, and again leads to a static probability density.

Ψ[x,y,3,1,1,3,α,β,t]/.{α->(π/4),β->(π/2)}

(E^((I π)/2))//ComplexExpand

Ψ2=(Sqrt[2]) (E^(-10 I (π^2) t)) (Sin[3 π x] Sin[π y]+ I Sin[π x] Sin[3 π y])

Only a single time factor remains, this time for the trivial reason

(3^2)+(1^2)=(1^2)+(3^2)=10

and as before, the probability density is static:

Ψ2star=(Sqrt[2]) (E^(+10 I (π^2) t)) (Sin[3 π x] Sin[π y]- I Sin[π x] Sin[3 π y])

Ψ2star*Ψ2//FullSimplify

P2[x_,y_]:=2 (((Sin[3 π x])^2) ((Sin[π y])^2)+((Sin[π x])^2) ((Sin[3 π y])^2))

Plot3D[P2[x,y],{x,0,1},{y,0,1}, Ticks->False]

Illustrative Case 3:

Ψ[x,y,3,1,1,2,α,β,t]/.{α->(π/4),β->(π/2)}

In this case

(3^2)+(1^2)=10
(1^1)+(2^2)=5

are distinct and so we can expect the probability density to move (oscillate).

Ψ3=(Sqrt[2]) ((E^(-10 I (π^2) t)) Sin[3 π x] Sin[π y]+ I (E^(-5 I (π^2) t))Sin[π x] Sin[2 π y])

Ψ3star=(Sqrt[2]) ((E^(+10 I (π^2) t)) Sin[3 π x] Sin[π y]- I (E^(+5 I (π^2) t))Sin[π x] Sin[2 π y])

Ψ3star*Ψ3//FullSimplify

P3[x_,y_,t_]:=2 (-4 (1+2 Cos[2 π x]) Cos[π y] Sin[5 (π^2) t] ((Sin[π x])^2) ((Sin[π y])^2)+((Sin[3 π x])^2) ((Sin[π y])^2)+((Sin[π x])^2) ((Sin[2 π y])^2))

This function is seen to oscillate with angular frequency 5(π^2), as shown below:

Manipulate[Plot3D[P3[x,y,((2π)/(50(π^2)))k],{x,0,1},{y,0,1}, Ticks->False, PlotRange->{0,6}], {k,0,10,1}]

The preceding examples (which I have borrowed from Mahoney, Paganin & Faulkner) will now be used  to expose the main points of this essay.

From probability current to quantum vorticity

Let us, for the purposes of the present argument, reinstate the physical parameters m and ℏ. Let  Ψ[x, y, t] be written in polar form

Ψ=R Exp[(I/ℏ)S];
Ψstar=R Exp[-(I/ℏ)S];

with

R=r[x,y,t];
S=s[x,y,t];

Supposing Ψ[x, y, t] to be a solution of the free-particle Schrödinger equation, we construct

Schrödinger=Exp[-(I/ℏ)S](-((ℏ^2)/(2m))(D[Ψ,{x,2}]+D[Ψ,{y,2}])-I ℏ D[Ψ,t])//Expand

from which we extract and massage the real and imaginary parts (which must vanish separately):

((CompConj[Schrödinger]+Schrödinger)/(2  r[x,y,t]))//Expand
((CompConj[Schrödinger]-Schrödinger)/(2 ℏ I))//Expand

The first equation can be written

Derivative[0,0,1][s][x,y,t]+(((Derivative[0,1,0][s][x,y,t])^2)/(2 m))+(((Derivative[1,0,0][s][x,y,t])^2)/(2 m))=(((ℏ^2)( Derivative[2,0,0][r][x,y,t]+Derivative[0,2,0][r][x,y,t]))/(2 m r[x,y,t]))

where in classical theory (ℏ ↓ 0) the expression on the left would vanish by appeal to Hamilton-Jacobi theory. The second equation finds employment at the end of the next paragraph.

The probability density is given now by

Ψstar *Ψ

and the probability current by

(Subscript[J, x])=(ℏ/(2m I))(Ψstar D[Ψ,x]-Ψ D[Ψstar,x])//Simplify
(Subscript[J, y])=(ℏ/(2m I))(Ψstar D[Ψ,y]-Ψ D[Ψstar,y])//Simplify

which can be expressed

probability current=(probability density)*((Grad[s])/m)

The continuity equation that provides a local account of probability conservation lends interest to a construction

D[((r[x,y,t])^2),t]+D[((((r[x,y,t])^2) Derivative[1,0,0][s][x,y,t])/m),x]+D[((((r[x,y,t])^2) Derivative[0,1,0][s][x,y,t])/m),y]//Expand

that after only a little massaging leads to a construction that has already been shown to vanish:

((%)/(r[x,y,t]))//Expand

%==((CompConj[Schrödinger]-Schrödinger)/( ℏ I))//Expand

The vector-valued construction

((probability current)/(probability density))

has the physical dimension of a "velocity." It becomes natural in this light to write

quantum mechanical velocity field ≡ ((probability current)/(probability density))=(1/m)Grad[S]

The velocity field

(({{(Subscript[V, x])[x,y,t]}, {(Subscript[V, y])[x,y,t]}}))=(((({{(Subscript[J, x])[x,y,t]}, {(Subscript[J, y])[x,y,t]}})))/(P[x,y,t]))=(1/m)Grad[S]

exists except where the probability density P[x,y,z,t] vanishes (which is to say: where Ψ[x,y,z,t] vanishes), and might exist even at such points if  J[x,y,z,t] happened also to vanish in just the right way. At those points the phase must be undifferentiable (discontinuous).

In fluid dynamics (see, for example, §1.2—"Rotation and Vorticity"—pages 19-31 in A. J. Chorin & J. E. Marsden, A Mathematical Introduction to Fluid Mechanics (2nd edition, 1990)) the vector-valued field produced by constructing the curl of a velocity field is called the "vorticity," denoated Ω:

Ω=Curl[V]

For the "quantum mechanical probability fluid" the vorticity vanishes

Ω=(1/m)Curl[ Grad[S]]=0

except at singular points (points where phase discontinuities cause  Grad[S]  to be undefined).

Concrete implications of the general theory

REMARK: We are working in two dimensions, so trhe probability current and Grad[s] are 2-vectors, while the vorticity Ω is effectively a scalar (z-component of a 3-vector).

Needs["VectorFieldPlots`"]

Application to Illustrative Case 1

This was the static case illustrated below:

The probability density vanishes at the center of the box

P1[(1/2),(1/2)]

where (therefore) the vorticity becomes singular. 

We construct and plot the probability current, which is seen to circulate (steadily) around the singularity:

-I(Ψ1star D[Ψ1,x]-Ψ1 D[Ψ1star,x])//ComplexExpand
-I(Ψ1star D[Ψ1,y]-Ψ1 D[Ψ1star,y])//ComplexExpand

J1X[x_,y_]:=-(64/9) π Cos[2 π x] Sin[π x] Sin[π y] Sin[2 π y]+(32/9) π Cos[π x] Sin[2 π x] Sin[π y] Sin[2 π y]
J1Y[x_,y_]:=(64/9) π Cos[2 π y] Sin[π x] Sin[2 π x] Sin[π y]-(32/9) π Cos[π y] Sin[π x] Sin[2 π x] Sin[2 π y]

VectorFieldPlot[{J1X[x,y],J1Y[x,y]},{x,0,1},{y,0,1},PlotLabel->"Steady Probability Current: Case 1", LabelStyle->Directive[Bold, Red]]

We construct and plot the associated velocity field (after modification of the denominator to keep Mathematica from complaining)

V1X[x_,y_]:=((J1X[x,y])/(P1[x,y]+0.000001))
V1Y[x_,y_]:=((J1Y[x,y])/(P1[x,y]+0.000001))

VectorFieldPlot[{V1X[x,y],V1Y[x,y]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field: Case 1", LabelStyle->Directive[Bold, Red]]

which places us in position to construct and plot the vorticity:

D[V1Y[x, y], x] - D[V1X[x, y], y]//Simplify//Chop

Ω1[x_,y_]:=(-0.00028073541407543063` ((Cos[π y])^2) ((Sin[π x])^3) Sin[π y]-0.00028073541407543063` ((Cos[π x])^2) Sin[π x] ((Sin[π y])^3)+0.00028073541407543063` ((Sin[π x])^3) ((Sin[π y])^3))/(((1.`*^-6+7.111111111111111` ((Cos[π x])^2) ((Sin[π x])^2) ((Sin[π y])^2)+7.111111111111111` Cos[π x] Cos[π y] ((Sin[π x])^2) ((Sin[π y])^2)+8.88888888888889` ((Cos[π y])^2) ((Sin[π x])^2) ((Sin[π y])^2)))^2)

Plot3D[Ω1[x,y],{x,0,1},{y,0,1}, PlotRange->{0,10},PlotPoints->{100}, Ticks->False,PlotLabel->"Vorticity: Case 1", LabelStyle->Directive[Bold, Red]]

The figure indicates that the vorticity vanishes everywhere except at the singularity, where it acquires the character of a delta function.

Though the amplitude of the wavefunction is (in the present example) time-independent

Ψ1

(as are all things derived from it: probability density, probability current, vorticity), its phase is time-dependent—with

Period=(2/(5π))

and has its own interesting story to tell:

Plot3D[Arg[ComplexExpand[WaveFunction1[x,y,0]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Initial Phase",Bold,Red]]

Plot3D[Arg[ComplexExpand[WaveFunction1[x,y,(Period/4)]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time (Period/4)",Bold,Red]]

Plot3D[Arg[ComplexExpand[WaveFunction1[x,y,(Period/2)]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time (Period/2)",Bold,Red]]

The complex phase exhibits a branch cut which ends at the singularity, and "twirls" in the sense indicated.

Application to Illustrative Case 2

This was again a static case, illustrated below:

The probability density now vanishes at four distinct points

P2[(1/3),(1/3)]
P2[(2/3),(1/3)]
P2[(1/3),(2/3)]
P2[(2/3),(2/3)]

We proceed exactly as before:

-I(Ψ2star D[Ψ2,x]-Ψ2 D[Ψ2star,x])//ComplexExpand
-I(Ψ2star D[Ψ2,y]-Ψ2 D[Ψ2star,y])//ComplexExpand

J2X[x_,y_]:=-12 π Cos[3 π x] Sin[π x] Sin[π y] Sin[3 π y]+4 π Cos[π x] Sin[3 π x] Sin[π y] Sin[3 π y]
J2Y[x_,y_]:=12 π Cos[3 π y] Sin[π x] Sin[3 π x] Sin[π y]-4 π Cos[π y] Sin[π x] Sin[3 π x] Sin[3 π y]

VectorFieldPlot[{J2X[x,y],J2Y[x,y]},{x,0,1},{y,0,1},PlotLabel->"Steady Probability Current: Case 2", LabelStyle->Directive[Bold, Red]]

V2X[x_,y_]:=((J2X[x,y])/(P2[x,y]+0.000001))
V2Y[x_,y_]:=((J2Y[x,y])/(P2[x,y]+0.000001))

VectorFieldPlot[{V2X[x,y],V2Y[x,y]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field: Case 2", LabelStyle->Directive[Bold, Red]]

D[V2Y[x, y], x] - D[V2X[x, y], y]//Simplify//Chop

Ω2[x_,y_]:=(0.0007106115168784338` Cos[3 π x] Cos[3 π y] Sin[π x] Sin[π y]-0.00007895683520871486` Cos[π x] Cos[π y] Sin[3 π x] Sin[3 π y])/(((1.`*^-6+2.` ((Sin[3 π x])^2) ((Sin[π y])^2)+2.` ((Sin[π x])^2) ((Sin[3 π y])^2)))^2)

Plot3D[Ω2[x,y],{x,.01,.99},{y,.01,.99}, PlotRange->{-10,10},PlotPoints->{100}, Ticks->False,PlotLabel->"Vorticity: Case 2", LabelStyle->Directive[Bold, Red]]

Ψ2

Period2=(1/(5π))

WaveFunction2[z_,y_,t_]:=(Sqrt[2]) (E^(-10 I (π^2) t)) (Sin[3 π x] Sin[π y]+I Sin[π x] Sin[3 π y])

Plot3D[Arg[ComplexExpand[WaveFunction2[x,y,0]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Initial Phase",Bold,Red]]

Plot3D[Arg[ComplexExpand[WaveFunction2[x,y,(Period2/4)]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time (Period/4)",Bold,Red]]

Plot3D[Arg[ComplexExpand[WaveFunction2[x,y,(Period2/2)]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time (Period/2)",Bold,Red]]
Plot3D[Arg[ComplexExpand[WaveFunction2[x,y,3(Period2/4)]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time 3(Period/4)",Bold,Red]]

Application to Illustrative Case 3

This case gains interest from the fact that the probability density

and associated current and vorticity are now time-dependent. I guess the position of an initial zero

P3[(2/3),(1/2),0]

(does the probability density vanish on a curve? hard to tell just by looking), which we know must move around cyclically.

Proceeding as before…

-I(Ψ1star D[Ψ3,x]-Ψ1 D[Ψ3star,x])//ComplexExpand//Simplify
-I(Ψ1star D[Ψ3,y]-Ψ1 D[Ψ3star,y])//ComplexExpand//Simplify

J3X[x_,y_,t_]:=(2/3) (Sqrt[2]) π Sin[2 π x] ((Sin[π y])^2) (2+6 Cos[π (5 π t-y)]-6 Cos[π (5 π t-2 x-y)]+4 Cos[π (x-y)]-6 Cos[π (5 π t+2 x-y)]+2 Cos[2 π y]+6 Cos[π (5 π t+y)]-6 Cos[π (5 π t-2 x+y)]+4 Cos[π (x+y)]-6 Cos[π (5 π t+2 x+y)]-6 Sin[π (5 π t-3 x)]-6 Sin[π (5 π t+3 x)]+3 Sin[π (5 π t-y)]-3 Sin[π (5 π t-2 x-y)]-3 Sin[π (5 π t+2 x-y)]+3 Sin[π (5 π t+y)]-3 Sin[π (5 π t-2 x+y)]-3 Sin[π (5 π t+2 x+y)])

J3Y[x_,y_,t_]:=(4/3) (Sqrt[2]) π (-Cos[5 (π^2) t] Cos[π y] Sin[3 π x] (2 Cos[10 (π^2) t] Sin[π x] Sin[2 π y]+Sin[10 (π^2) t] (2 Sin[2 π x] Sin[π y]+Sin[π x] Sin[2 π y]))+2 ((Cos[5 (π^2) t])^2) (2 Cos[π x]+Cos[π y]) ((Sin[π x])^2) (-Sin[π y]+Sin[3 π y])+Sin[5 (π^2) t] (4 Cos[2 π y] Sin[5 (π^2) t] Sin[π x] Sin[2 π x] Sin[π y]+Cos[10 (π^2) t] Cos[π y] Sin[3 π x] (2 Sin[2 π x] Sin[π y]+Sin[π x] Sin[2 π y])+Sin[π x] (-2 Cos[π y] Sin[10 (π^2) t] Sin[3 π x] Sin[2 π y]+Sin[5 (π^2) t] Sin[π x] Sin[4 π y])))

Evidently the current configuration oscillates with

Period3=(2/(5π))

VectorFieldPlot[{J3X[x,y,0],J3Y[x,y,0]},{x,0,1},{y,0,1},PlotLabel->"Initial Probability Current: Case 3", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{J3X[x,y,(Period3/16)],J3Y[x,y,(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Probability Current at Time (Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{J3X[x,y,2(Period3/16)],J3Y[x,y,2(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Probability Current at Time 2(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{J3X[x,y,3(Period3/16)],J3Y[x,y,3(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Probability Current at Time 3(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{J3X[x,y,4(Period3/16)],J3Y[x,y,4(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Probability Current at Time 4(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{J3X[x,y,5(Period3/16)],J3Y[x,y,5(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Probability Current at Time 5(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{J3X[x,y,6(Period3/16)],J3Y[x,y,6(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Probability Current at Time 6(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{J3X[x,y,7(Period3/16)],J3Y[x,y,7(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Probability Current at Time 7(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{J3X[x,y,8(Period3/16)],J3Y[x,y,8(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Probability Current at Time 8(Period/16)", LabelStyle->Directive[Bold, Red]]

V3X[x_,y_,t_]:=((J3X[x,y,t])/(P3[x,y,t]+0.000001))
V3Y[x_,y_,t_]:=((J3Y[x,y,t])/(P3[x,y,t]+0.000001))

VectorFieldPlot[{V3X[x,y,0],V3Y[x,y,0]},{x,0,1},{y,0,1},PlotLabel->"Initial Velocity Field: Case 3", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{V3X[x,y,(Period3/16)],V3Y[x,y,(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field at Time (Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{V3X[x,y,2(Period3/16)],V3Y[x,y,2(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field at Time 2(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{V3X[x,y,3(Period3/16)],V3Y[x,y,3(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field at Time 3(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{V3X[x,y,4(Period3/16)],V3Y[x,y,4(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field at Time 4(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{V3X[x,y,5(Period3/16)],V3Y[x,y,5(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field at Time 5(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{V3X[x,y,6(Period3/16)],V3Y[x,y,6(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field at Time 6(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{V3X[x,y,7(Period3/16)],V3Y[x,y,7(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field at Time 7(Period/16)", LabelStyle->Directive[Bold, Red]]

VectorFieldPlot[{V3X[x,y,8(Period3/16)],V3Y[x,y,8(Period3/16)]},{x,0,1},{y,0,1},PlotLabel->"Velocity Field at Time 8(Period/16)", LabelStyle->Directive[Bold, Red]]

We look finally to the motion of the complex phase:

Ψ3

WaveFunction3[z_,y_,t_]:=(Sqrt[2]) ((E^(-10 I (π^2) t)) Sin[3 π x] Sin[π y]+I (E^(-5 I (π^2) t)) Sin[π x] Sin[2 π y])

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunction3[x,y,0]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Initial Phase",Bold,Red]]

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunction3[x,y,(Period3/16)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time (Period/16)",Bold,Red]]
Plot3D[Evaluate[Arg[ComplexExpand[WaveFunction3[x,y,2(Period3/16)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time 2(Period/16)",Bold,Red]]
Plot3D[Evaluate[Arg[ComplexExpand[WaveFunction3[x,y,3(Period3/16)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time 3(Period/16)",Bold,Red]]

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunction3[x,y,4(Period3/16)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time 4(Period/16)",Bold,Red]]
Plot3D[Evaluate[Arg[ComplexExpand[WaveFunction3[x,y,5(Period3/16)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time 5(Period/16)",Bold,Red]]
Plot3D[Evaluate[Arg[ComplexExpand[WaveFunction3[x,y,6(Period3/16)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time 6(Period/16)",Bold,Red]]
Plot3D[Evaluate[Arg[ComplexExpand[WaveFunction3[x,y,7(Period3/16)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time 7(Period/16)",Bold,Red]]
Plot3D[Evaluate[Arg[ComplexExpand[WaveFunction3[x,y,8(Period3/16)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}},PlotLabel->Style["Phase at time 8(Period/16)",Bold,Red]]

Application to the case of a simple quantum Lissajous

We have demonstrated that the wavefunction

Φ=(1.1343346991218932` (E^(-I t)) Sin[π x]+0.8445618925620814` (E^(-4 I t)) Sin[2 π x]) (-1.144048972734524` (E^(-9 I t)) Sin[3 π y]-0.8313554883352129` (E^(-16 I t)) Sin[4 π y])

Φstar=CompConj[Φ]

leads to a probability density that moves with hardly more definition than a cloud

causes expected value of {x,y} to trace a nice crisp Lissajous figure of the form

When we look to the probability current, we encounter expressions that are in general very complicated

-I(Φstar D[Φ,x]-Φ D[Φstar,x])//ComplexExpand//Chop
-I(Φstar D[Φ,y]-Φ D[Φstar,y])//ComplexExpand//Chop

JLissX[x_,y_,t_]:=15.756936833748924` Cos[4 t] ((Cos[9 t])^2) Cos[2 π x] Sin[t] Sin[π x] ((Sin[3 π y])^2)-15.756936833748924` Cos[t] ((Cos[9 t])^2) Cos[2 π x] Sin[4 t] Sin[π x] ((Sin[3 π y])^2)+15.756936833748924` Cos[4 t] Cos[2 π x] Sin[t] ((Sin[9 t])^2) Sin[π x] ((Sin[3 π y])^2)-15.756936833748924` Cos[t] Cos[2 π x] Sin[4 t] ((Sin[9 t])^2) Sin[π x] ((Sin[3 π y])^2)-7.878468416874462` Cos[4 t] ((Cos[9 t])^2) Cos[π x] Sin[t] Sin[2 π x] ((Sin[3 π y])^2)+7.878468416874462` Cos[t] ((Cos[9 t])^2) Cos[π x] Sin[4 t] Sin[2 π x] ((Sin[3 π y])^2)-7.878468416874462` Cos[4 t] Cos[π x] Sin[t] ((Sin[9 t])^2) Sin[2 π x] ((Sin[3 π y])^2)+7.878468416874462` Cos[t] Cos[π x] Sin[4 t] ((Sin[9 t])^2) Sin[2 π x] ((Sin[3 π y])^2)+22.900446096774214` Cos[4 t] Cos[9 t] Cos[16 t] Cos[2 π x] Sin[t] Sin[π x] Sin[3 π y] Sin[4 π y]-22.900446096774214` Cos[t] Cos[9 t] Cos[16 t] Cos[2 π x] Sin[4 t] Sin[π x] Sin[3 π y] Sin[4 π y]+22.900446096774214` Cos[4 t] Cos[2 π x] Sin[t] Sin[9 t] Sin[16 t] Sin[π x] Sin[3 π y] Sin[4 π y]-22.900446096774214` Cos[t] Cos[2 π x] Sin[4 t] Sin[9 t] Sin[16 t] Sin[π x] Sin[3 π y] Sin[4 π y]-11.450223048387107` Cos[4 t] Cos[9 t] Cos[16 t] Cos[π x] Sin[t] Sin[2 π x] Sin[3 π y] Sin[4 π y]+11.450223048387107` Cos[t] Cos[9 t] Cos[16 t] Cos[π x] Sin[4 t] Sin[2 π x] Sin[3 π y] Sin[4 π y]-11.450223048387107` Cos[4 t] Cos[π x] Sin[t] Sin[9 t] Sin[16 t] Sin[2 π x] Sin[3 π y] Sin[4 π y]+11.450223048387107` Cos[t] Cos[π x] Sin[4 t] Sin[9 t] Sin[16 t] Sin[2 π x] Sin[3 π y] Sin[4 π y]+8.320627875908158` Cos[4 t] ((Cos[16 t])^2) Cos[2 π x] Sin[t] Sin[π x] ((Sin[4 π y])^2)-8.320627875908158` Cos[t] ((Cos[16 t])^2) Cos[2 π x] Sin[4 t] Sin[π x] ((Sin[4 π y])^2)+8.320627875908158` Cos[4 t] Cos[2 π x] Sin[t] ((Sin[16 t])^2) Sin[π x] ((Sin[4 π y])^2)-8.320627875908158` Cos[t] Cos[2 π x] Sin[4 t] ((Sin[16 t])^2) Sin[π x] ((Sin[4 π y])^2)-4.160313937954079` Cos[4 t] ((Cos[16 t])^2) Cos[π x] Sin[t] Sin[2 π x] ((Sin[4 π y])^2)+4.160313937954079` Cos[t] ((Cos[16 t])^2) Cos[π x] Sin[4 t] Sin[2 π x] ((Sin[4 π y])^2)-4.160313937954079` Cos[4 t] Cos[π x] Sin[t] ((Sin[16 t])^2) Sin[2 π x] ((Sin[4 π y])^2)+4.160313937954079` Cos[t] Cos[π x] Sin[4 t] ((Sin[16 t])^2) Sin[2 π x] ((Sin[4 π y])^2)

JLissY[x_,y_,t_]:=30.7576873426503` ((Cos[t])^2) Cos[16 t] Cos[4 π y] Sin[9 t] ((Sin[π x])^2) Sin[3 π y]+30.7576873426503` Cos[16 t] Cos[4 π y] ((Sin[t])^2) Sin[9 t] ((Sin[π x])^2) Sin[3 π y]-30.7576873426503` ((Cos[t])^2) Cos[9 t] Cos[4 π y] Sin[16 t] ((Sin[π x])^2) Sin[3 π y]-30.7576873426503` Cos[9 t] Cos[4 π y] ((Sin[t])^2) Sin[16 t] ((Sin[π x])^2) Sin[3 π y]+45.80089219354843` Cos[t] Cos[4 t] Cos[16 t] Cos[4 π y] Sin[9 t] Sin[π x] Sin[2 π x] Sin[3 π y]+45.80089219354843` Cos[16 t] Cos[4 π y] Sin[t] Sin[4 t] Sin[9 t] Sin[π x] Sin[2 π x] Sin[3 π y]-45.80089219354843` Cos[t] Cos[4 t] Cos[9 t] Cos[4 π y] Sin[16 t] Sin[π x] Sin[2 π x] Sin[3 π y]-45.80089219354843` Cos[9 t] Cos[4 π y] Sin[t] Sin[4 t] Sin[16 t] Sin[π x] Sin[2 π x] Sin[3 π y]+17.050385667457427` ((Cos[4 t])^2) Cos[16 t] Cos[4 π y] Sin[9 t] ((Sin[2 π x])^2) Sin[3 π y]+17.050385667457427` Cos[16 t] Cos[4 π y] ((Sin[4 t])^2) Sin[9 t] ((Sin[2 π x])^2) Sin[3 π y]-17.050385667457427` ((Cos[4 t])^2) Cos[9 t] Cos[4 π y] Sin[16 t] ((Sin[2 π x])^2) Sin[3 π y]-17.050385667457427` Cos[9 t] Cos[4 π y] ((Sin[4 t])^2) Sin[16 t] ((Sin[2 π x])^2) Sin[3 π y]-23.06826550698773` ((Cos[t])^2) Cos[16 t] Cos[3 π y] Sin[9 t] ((Sin[π x])^2) Sin[4 π y]-23.06826550698773` Cos[16 t] Cos[3 π y] ((Sin[t])^2) Sin[9 t] ((Sin[π x])^2) Sin[4 π y]+23.06826550698773` ((Cos[t])^2) Cos[9 t] Cos[3 π y] Sin[16 t] ((Sin[π x])^2) Sin[4 π y]+23.06826550698773` Cos[9 t] Cos[3 π y] ((Sin[t])^2) Sin[16 t] ((Sin[π x])^2) Sin[4 π y]-34.35066914516132` Cos[t] Cos[4 t] Cos[16 t] Cos[3 π y] Sin[9 t] Sin[π x] Sin[2 π x] Sin[4 π y]-34.35066914516132` Cos[16 t] Cos[3 π y] Sin[t] Sin[4 t] Sin[9 t] Sin[π x] Sin[2 π x] Sin[4 π y]+34.35066914516132` Cos[t] Cos[4 t] Cos[9 t] Cos[3 π y] Sin[16 t] Sin[π x] Sin[2 π x] Sin[4 π y]+34.35066914516132` Cos[9 t] Cos[3 π y] Sin[t] Sin[4 t] Sin[16 t] Sin[π x] Sin[2 π x] Sin[4 π y]-12.78778925059307` ((Cos[4 t])^2) Cos[16 t] Cos[3 π y] Sin[9 t] ((Sin[2 π x])^2) Sin[4 π y]-12.78778925059307` Cos[16 t] Cos[3 π y] ((Sin[4 t])^2) Sin[9 t] ((Sin[2 π x])^2) Sin[4 π y]+12.78778925059307` ((Cos[4 t])^2) Cos[9 t] Cos[3 π y] Sin[16 t] ((Sin[2 π x])^2) Sin[4 π y]+12.78778925059307` Cos[9 t] Cos[3 π y] ((Sin[4 t])^2) Sin[16 t] ((Sin[2 π x])^2) Sin[4 π y]

—though at certain times they do assume greatly simpified form

JLissX[x,y,0]
JLissY[x,y,0]

JLissX[x,y,π]
JLissY[x,y,π]

JLissX[x,y,1]
JLissY[x,y,1]

—and which give rise to not-very-informative vector field plots:

Table[{k,VectorFieldPlot[Evaluate[{JLissX[x,y,(k/5)+.0001],JLissY[x,y,(k/5)+.0001]}],{x,0,1},{y,0,1}]},{k,0,25}]//MatrixForm

When we look finally to the phase of the wavefunction,

Φ

ComplexExpand[Φ]

WaveFunctionΦ[x_,y_,t_]:=(1.1343346991218932` (E^(-I t)) Sin[π x]+0.8445618925620814` (E^(-4 I t)) Sin[2 π x]) (-1.144048972734524` (E^(-9 I t)) Sin[3 π y]-0.8313554883352129` (E^(-16 I t)) Sin[4 π y])

we are led to figures

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunctionΦ[x,y,(1/5)k]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},  Ticks->{{},{},{-π,π}}]

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunctionΦ[x,y,(0/5)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},PlotPoints->20 , Ticks->{{},{},{-π,π}}]

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunctionΦ[x,y,(1/5)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},PlotPoints->20 , Ticks->{{},{},{-π,π}}]

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunctionΦ[x,y,(2/5)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},PlotPoints->20 , Ticks->{{},{},{-π,π}}]

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunctionΦ[x,y,(3/5)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},PlotPoints->20 , Ticks->{{},{},{-π,π}}]

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunctionΦ[x,y,(4/5)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},PlotPoints->20 , Ticks->{{},{},{-π,π}}]

Plot3D[Evaluate[Arg[ComplexExpand[WaveFunctionΦ[x,y,(5/5)]]]],{x,0,1},{y,0,1}, PlotRange->{-π,π},PlotPoints->20 , Ticks->{{},{},{-π,π}}]

that are complicated at every instant, and that move in complicated ways.

Appliction to Chen's quantum Lissajous state

We had

Ψ[x_,y_,m_,n_]:=2 Sin[m π x]Sin[n π y]
Chen[x_,y_,n_, p_, q_,ϕ_]:=(1/(2^(n/2)))(∑_(k=0)^n)(Sqrt[Binomial[n,k]])(E^(I k ϕ))Ψ[p x,q y,k+1,n-k+1]

which led to real-valued intial state like this

Chen[x,y,10, 7,3,0]

ChenΨ[x_,y_]:=(1/32) (2 Sin[77 π x] Sin[3 π y]+2 (Sqrt[10]) Sin[70 π x] Sin[6 π y]+6 (Sqrt[5]) Sin[63 π x] Sin[9 π y]+4 (Sqrt[30]) Sin[56 π x] Sin[12 π y]+2 (Sqrt[210]) Sin[49 π x] Sin[15 π y]+12 (Sqrt[7]) Sin[42 π x] Sin[18 π y]+2 (Sqrt[210]) Sin[35 π x] Sin[21 π y]+4 (Sqrt[30]) Sin[28 π x] Sin[24 π y]+6 (Sqrt[5]) Sin[21 π x] Sin[27 π y]+2 (Sqrt[10]) Sin[14 π x] Sin[30 π y]+2 Sin[7 π x] Sin[33 π y])

and which led to probability density plots like this:

When we look to the definition of probability current

JX=-i(Ψstar D[Ψ,x]-Ψ D[Ψstar,x])//ComplexExpand
JY=-i(Ψstar D[Ψ,y]-Ψ D[Ψstar,y])//ComplexExpand

we see that to avoid the uninteresting statements

JX =0
JY =0

the state Ψ must be complex. The Chen state to which we look previously was, however, real, so we look now to the Chen state

Chen[x,y,10, 7,3,(π/3)]//ComplexExpand

ChenΨ=-(1/32) Sin[77 π x] Sin[3 π y]-(1/8) (Sqrt[(5/2)]) Sin[70 π x] Sin[6 π y]-(3/32) (Sqrt[5]) Sin[63 π x] Sin[9 π y]+(1/8) (Sqrt[(15/2)]) Sin[56 π x] Sin[12 π y]+(1/8) (Sqrt[(105/2)]) Sin[49 π x] Sin[15 π y]+(3/16) (Sqrt[7]) Sin[42 π x] Sin[18 π y]-(1/16) (Sqrt[(105/2)]) Sin[35 π x] Sin[21 π y]-(1/4) (Sqrt[(15/2)]) Sin[28 π x] Sin[24 π y]-(3/32) (Sqrt[5]) Sin[21 π x] Sin[27 π y]+(1/16) (Sqrt[(5/2)]) Sin[14 π x] Sin[30 π y]+I (-(1/32) (Sqrt[3]) Sin[77 π x] Sin[3 π y]+(3/32) (Sqrt[15]) Sin[63 π x] Sin[9 π y]+(3/8) (Sqrt[(5/2)]) Sin[56 π x] Sin[12 π y]-(3/16) (Sqrt[21]) Sin[42 π x] Sin[18 π y]-(3/16) (Sqrt[(35/2)]) Sin[35 π x] Sin[21 π y]+(3/32) (Sqrt[15]) Sin[21 π x] Sin[27 π y]+(1/16) (Sqrt[(15/2)]) Sin[14 π x] Sin[30 π y])+(1/16) Sin[7 π x] Sin[33 π y];

RealChenΨ=-(1/32) Sin[77 π x] Sin[3 π y]-(1/8) (Sqrt[(5/2)]) Sin[70 π x] Sin[6 π y]-(3/32) (Sqrt[5]) Sin[63 π x] Sin[9 π y]+(1/8) (Sqrt[(15/2)]) Sin[56 π x] Sin[12 π y]+(1/8) (Sqrt[(105/2)]) Sin[49 π x] Sin[15 π y]+(3/16) (Sqrt[7]) Sin[42 π x] Sin[18 π y]-(1/16) (Sqrt[(105/2)]) Sin[35 π x] Sin[21 π y]-(1/4) (Sqrt[(15/2)]) Sin[28 π x] Sin[24 π y]-(3/32) (Sqrt[5]) Sin[21 π x] Sin[27 π y]+(1/16) (Sqrt[(5/2)]) Sin[14 π x] Sin[30 π y]+(1/16) Sin[7 π x] Sin[33 π y];
ImagChenΨ=-(1/32) (Sqrt[3]) Sin[77 π x] Sin[3 π y]+(3/32) (Sqrt[15]) Sin[63 π x] Sin[9 π y]+(3/8) (Sqrt[(5/2)]) Sin[56 π x] Sin[12 π y]-(3/16) (Sqrt[21]) Sin[42 π x] Sin[18 π y]-(3/16) (Sqrt[(35/2)]) Sin[35 π x] Sin[21 π y]+(3/32) (Sqrt[15]) Sin[21 π x] Sin[27 π y]+(1/16) (Sqrt[(15/2)]) Sin[14 π x] Sin[30 π y];

DensityPlot[Evaluate[((RealChenΨ)^2)+((ImagChenΨ)^2)],{x,0,1},{y,0,1}, PlotPoints->100]

DensityPlot[Evaluate[((RealChenΨ)^2)+((ImagChenΨ)^2)],{x,0,(4/7)},{y,0,(2/3)}, PlotPoints->100]

ChenΨstar=RealChenΨ-I ImagChenΨ;

JXChen=-I(ChenΨstar D[ChenΨ,x]-ChenΨ D[ChenΨstar,x])//ComplexExpand;
JYChen=-I(ChenΨstar D[ChenΨ,y]-ChenΨ D[ChenΨstar,y])//ComplexExpand;

Length[JXChen]
Length[JYChen]

VectorFieldPlot[{JXChen,JYChen},{x,0,(4/7)},{y,0,(2/3)}, PlotPoints->40]

VectorFieldPlot[{JXChen,JYChen},{x,0,(1/7)},{y,0,(1/3)}, PlotPoints->9]

VectorFieldPlot[{JXChen,JYChen},{x,0,(2/7),(1/28)},{y,0,(1/3), (1/30)}]